\chapter{Interpolation and Smoothing Algorithms}
\label{ch:interpolation_smoothing}

\section{Introduction}

Interpolation and smoothing algorithms form a critical component of the GSI (Gridpoint Statistical Interpolation) system, enabling the accurate transfer of information between different spatial representations and the application of spatial filters for numerical stability and quality control. These algorithms must handle diverse data types including scalar fields (temperature, pressure, humidity), vector fields (wind components), and specialized variables (ozone, aerosols) while maintaining physical consistency and numerical accuracy.

The challenge lies in developing interpolation schemes that preserve important physical properties such as mass conservation, monotonicity, and bounded variation while providing the computational efficiency required for operational data assimilation systems. This chapter examines the comprehensive suite of interpolation and smoothing routines implemented in GSI, covering both traditional methods and advanced algorithms tailored for atmospheric applications.

\section{Fundamental Interpolation Methods}

\subsection{Bilinear Analysis Grid Interpolation}

The core analysis grid interpolation routine (\texttt{intrp2a.f90}) implements sophisticated bilinear interpolation algorithms optimized for atmospheric data assimilation. Unlike simple bilinear interpolation, this implementation accounts for:

\begin{itemize}
    \item Spherical geometry effects
    \item Terrain-following coordinates
    \item Variable density grids
    \item Boundary condition handling
    \item Conservation properties
\end{itemize}

The basic bilinear interpolation formula is enhanced for atmospheric applications:

\begin{equation}
\phi(x, y) = \sum_{i,j} w_{i,j}(x, y) \phi_{i,j}
\end{equation}

where the weights $w_{i,j}$ are computed considering the spherical Earth geometry:

\begin{align}
w_{i,j} &= \frac{A_{i,j}(x, y)}{A_{total}} \\
A_{i,j} &= R^2 \cos(\phi_{ij}) \Delta\lambda \Delta\phi
\end{align}

Here $R$ is Earth's radius, $\phi_{ij}$ is the latitude at grid point $(i,j)$, and $\Delta\lambda$, $\Delta\phi$ are the longitude and latitude increments.

\textbf{Monotonicity Preservation:}
To maintain monotonicity in critical atmospheric variables, the algorithm employs flux-corrected transport principles:

\begin{equation}
\phi_{new} = \phi_{old} + \alpha \cdot \min(\phi_{max} - \phi_{old}, \phi_{old} - \phi_{min}, \Delta\phi_{uncorrected})
\end{equation}

where $\alpha$ is a limiter coefficient and $\phi_{max}$, $\phi_{min}$ are local extrema bounds.

\subsection{Three-Dimensional Ozone Interpolation}

The three-dimensional ozone interpolation routine (\texttt{intrp3oz.f90}) addresses the unique challenges of ozone distribution interpolation in the stratosphere and troposphere. Ozone exhibits strong vertical gradients, seasonal variations, and chemical-dynamical coupling that require specialized treatment.

The algorithm employs a hybrid approach combining:

\textbf{Vertical Coordinate Transformation:}
\begin{equation}
z^* = \frac{\ln(p_{surface}) - \ln(p)}{\ln(p_{surface}) - \ln(p_{top})}
\end{equation}

This log-pressure coordinate provides better resolution in the stratosphere where ozone concentrations are highest.

\textbf{Spline-Based Vertical Interpolation:}
\begin{align}
O_3(z^*) &= \sum_{k=0}^{3} c_k (z^* - z^*_k)^k \\
c_k &= \text{coefficients from cubic spline fitting}
\end{align}

The cubic splines preserve the sharp gradients at the tropopause while providing smooth interpolation in regions with sparse observations.

\textbf{Horizontal Anisotropic Interpolation:}
Ozone fields exhibit strong latitudinal gradients, requiring anisotropic interpolation:

\begin{equation}
O_3(\lambda, \phi) = \sum_{i,j} w_{i,j}(\lambda, \phi) O_{3,i,j} \cdot f(\phi)
\end{equation}

where $f(\phi)$ is a latitude-dependent scaling factor accounting for the poleward ozone increase.

\subsection{Masked Field Interpolation}

The masked field interpolation routine (\texttt{intrp\_msk.f90}) handles interpolation in the presence of missing data, land-sea boundaries, and terrain obstacles. This is particularly important for surface and near-surface observations where topography and surface type create data voids.

The algorithm employs a multi-stage approach:

\textbf{Stage 1: Valid Neighbor Search}
\begin{equation}
N_{valid}(x,y) = \{(i,j) : |i-x| \leq R, |j-y| \leq R, mask_{i,j} = 1\}
\end{equation}

where $R$ is the search radius and $mask_{i,j}$ indicates valid data points.

\textbf{Stage 2: Distance-Weighted Interpolation}
\begin{equation}
\phi(x,y) = \frac{\sum_{(i,j) \in N_{valid}} w_{i,j} \phi_{i,j}}{\sum_{(i,j) \in N_{valid}} w_{i,j}}
\end{equation}

with inverse distance weighting:
\begin{equation}
w_{i,j} = \frac{1}{(d_{i,j}^2 + \epsilon)^p}
\end{equation}

where $d_{i,j}$ is the distance to grid point $(i,j)$, $\epsilon$ is a small regularization parameter, and $p$ controls the sharpness of the weighting (typically $p = 1$ or $p = 2$).

\textbf{Stage 3: Boundary Condition Enforcement}
Near boundaries, the algorithm applies:
\begin{equation}
\phi_{boundary} = \phi_{interior} + \nabla\phi \cdot \vec{n} \cdot d
\end{equation}

where $\vec{n}$ is the outward normal and $d$ is the distance to the boundary.

\section{Temporal Interpolation Methods}

\subsection{Time-Based Analysis Interpolation}

The temporal interpolation routine (\texttt{tintrp2a.f90}) addresses the temporal aspect of data assimilation, interpolating background fields and observations to the analysis time. This is crucial for 4DVar systems and for handling asynchronous observations.

The algorithm implements several temporal interpolation schemes:

\textbf{Linear Temporal Interpolation:}
\begin{equation}
\phi(t) = \phi(t_1) + \frac{t - t_1}{t_2 - t_1}[\phi(t_2) - \phi(t_1)]
\end{equation}

\textbf{Quadratic Temporal Interpolation:}
For better accuracy with rapidly evolving fields:
\begin{align}
\phi(t) &= \phi(t_1) + a(t - t_1) + b(t - t_1)^2 \\
a &= \frac{\phi(t_2) - \phi(t_0)}{t_2 - t_0} \\
b &= \frac{\phi(t_2) - 2\phi(t_1) + \phi(t_0)}{(t_2 - t_0)^2/4}
\end{align}

\textbf{Advective Temporal Interpolation:}
For tracer species and moisture fields:
\begin{equation}
\phi(x, t) = \phi(x - \vec{u} \Delta t, t_0)
\end{equation}

where $\vec{u}$ is the advection velocity field.

\subsection{Three-Dimensional Temporal Interpolation}

The three-dimensional temporal interpolation routine (\texttt{tintrp3.f90}) extends temporal interpolation to three-dimensional fields while maintaining vertical coordinate consistency and physical constraints.

The algorithm addresses:

\begin{itemize}
    \item Pressure coordinate variations with time
    \item Diabatic heating effects on vertical motion
    \item Boundary layer evolution
    \item Mass conservation in time
\end{itemize}

\textbf{Isentropic Temporal Interpolation:}
For upper-level fields, the algorithm uses potential temperature as a vertical coordinate:
\begin{equation}
\phi(\theta, t) = \phi(\theta, t_1) + \frac{t - t_1}{t_2 - t_1}[\phi(\theta, t_2) - \phi(\theta, t_1)]
\end{equation}

This preserves adiabatic evolution and reduces interpolation errors in stratified flows.

\textbf{Boundary Layer Temporal Interpolation:}
In the boundary layer, the algorithm accounts for diurnal cycles:
\begin{equation}
\phi(z, t) = \phi_{mean}(z) + A(z)\cos(\omega t + \phi_0(z))
\end{equation}

where $\omega = 2\pi/24$ hours and $A(z)$, $\phi_0(z)$ are height-dependent amplitude and phase functions.

\section{Spatial Smoothing Algorithms}

\subsection{Polar Coordinate Smoothing}

The polar coordinate smoothing routine (\texttt{smooth\_polcarf.f90}) addresses the unique computational challenges near the poles in global atmospheric models. Standard Cartesian smoothing operators become ill-conditioned near the poles due to grid convergence, requiring specialized formulations.

The algorithm transforms the smoothing operator to polar coordinates:

\textbf{Laplacian in Polar Coordinates:}
\begin{equation}
\nabla^2 \phi = \frac{1}{r^2}\left[\frac{\partial}{\partial r}\left(r^2 \frac{\partial \phi}{\partial r}\right) + \frac{1}{\sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta \frac{\partial \phi}{\partial \theta}\right) + \frac{1}{\sin^2\theta}\frac{\partial^2 \phi}{\partial \lambda^2}\right]
\end{equation}

Near the poles ($\theta \approx 0$ or $\theta \approx \pi$), the algorithm uses:

\textbf{Polar Cap Treatment:}
\begin{equation}
\nabla^2 \phi|_{pole} = \frac{2}{R^2}\left[\frac{1}{N}\sum_{k=1}^{N} \phi(\lambda_k, \theta_{near}) - \phi(pole)\right]
\end{equation}

where the sum is over $N$ equally spaced longitude points at the nearest latitude circle.

\textbf{Smoothing Filter Application:}
The smoothing operator is applied iteratively:
\begin{equation}
\phi^{n+1} = \phi^n + \alpha \nabla^2 \phi^n
\end{equation}

with stability condition $\alpha < \Delta t_{max}$ where:
\begin{equation}
\Delta t_{max} = \frac{\min(\Delta x^2, \Delta y^2)}{4D}
\end{equation}

and $D$ is the diffusion coefficient.

\subsection{WRF-Specific Wind Smoothing}

The WRF wind smoothing routine (\texttt{smoothwwrf.f90}) applies specialized smoothing to wind components in Weather Research and Forecasting (WRF) model grids. The C-grid staggering of WRF requires careful treatment to maintain:

\begin{itemize}
    \item Geostrophic balance
    \item Mass continuity
    \item Boundary layer structure
    \item Topographic flow features
\end{itemize}

\textbf{Vector Field Smoothing:}
The algorithm smooths the vorticity and divergence components separately:

\begin{align}
\zeta &= \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \\
\delta &= \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}
\end{align}

After applying smoothing operators $S$ to each component:
\begin{align}
\zeta_{smooth} &= S[\zeta] \\
\delta_{smooth} &= S[\delta]
\end{align}

The wind components are reconstructed using a stream function $\psi$ and velocity potential $\chi$:

\begin{align}
u &= -\frac{\partial \psi}{\partial y} + \frac{\partial \chi}{\partial x} \\
v &= \frac{\partial \psi}{\partial x} + \frac{\partial \chi}{\partial y}
\end{align}

where:
\begin{align}
\nabla^2 \psi &= \zeta_{smooth} \\
\nabla^2 \chi &= \delta_{smooth}
\end{align}

\textbf{Geostrophic Balance Preservation:}
To maintain geostrophic balance, the algorithm applies the constraint:
\begin{equation}
\vec{u}_g = \frac{1}{f}\vec{k} \times \nabla \Phi
\end{equation}

where $f$ is the Coriolis parameter and $\Phi$ is the geopotential.

\subsection{Generalized Spatial Smoothing}

The generalized spatial smoothing routine (\texttt{smoothzrf.f90}) provides a flexible framework for applying various smoothing operators to three-dimensional atmospheric fields. The routine supports:

\begin{itemize}
    \item Isotropic and anisotropic smoothing
    \item Level-dependent smoothing parameters
    \item Variable-specific smoothing coefficients
    \item Boundary-aware smoothing near terrain
\end{itemize}

\textbf{Anisotropic Smoothing Tensor:}
\begin{equation}
S[\phi] = \nabla \cdot (D \nabla \phi)
\end{equation}

where $D$ is the anisotropic diffusion tensor:
\begin{equation}
D = \begin{pmatrix}
D_{xx} & D_{xy} \\
D_{xy} & D_{yy}
\end{pmatrix}
\end{equation}

For atmospheric applications, the tensor components are typically:
\begin{align}
D_{xx} &= D_0 \cos^2(\alpha) + D_1 \sin^2(\alpha) \\
D_{yy} &= D_0 \sin^2(\alpha) + D_1 \cos^2(\alpha) \\
D_{xy} &= (D_0 - D_1) \sin(\alpha) \cos(\alpha)
\end{align}

where $\alpha$ is the angle of principal diffusion axis and $D_0$, $D_1$ are the principal diffusion coefficients.

\section{Specialized Processing Algorithms}

\subsection{Derivative Computation}

The derivative computation routines (\texttt{get\_derivatives.f90}, \texttt{get\_derivatives2.f90}) calculate spatial derivatives of atmospheric fields using finite difference approximations optimized for atmospheric applications.

\textbf{First-Order Derivatives:}
Standard second-order accurate finite differences:
\begin{align}
\frac{\partial \phi}{\partial x} &\approx \frac{\phi_{i+1,j} - \phi_{i-1,j}}{2\Delta x} \\
\frac{\partial \phi}{\partial y} &\approx \frac{\phi_{i,j+1} - \phi_{i,j-1}}{2\Delta y}
\end{align}

\textbf{Higher-Order Derivatives:}
Fourth-order accurate approximations for smooth fields:
\begin{align}
\frac{\partial \phi}{\partial x} &\approx \frac{-\phi_{i+2,j} + 8\phi_{i+1,j} - 8\phi_{i-1,j} + \phi_{i-2,j}}{12\Delta x} \\
\frac{\partial^2 \phi}{\partial x^2} &\approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\Delta x^2}
\end{align}

\textbf{Spherical Coordinate Derivatives:}
On the sphere, derivatives account for the metric terms:
\begin{align}
\frac{\partial \phi}{\partial \lambda} &\approx \frac{\phi(\lambda + \Delta\lambda) - \phi(\lambda - \Delta\lambda)}{2R\cos\theta \Delta\lambda} \\
\frac{\partial \phi}{\partial \theta} &\approx \frac{\phi(\theta + \Delta\theta) - \phi(\theta - \Delta\theta)}{2R\Delta\theta}
\end{align}

\textbf{Mixed Derivatives:}
Cross derivatives for computing vorticity and shear:
\begin{equation}
\frac{\partial^2 \phi}{\partial x \partial y} \approx \frac{\phi_{i+1,j+1} - \phi_{i+1,j-1} - \phi_{i-1,j+1} + \phi_{i-1,j-1}}{4\Delta x \Delta y}
\end{equation}

\subsection{Pressure Field Interpolation}

The pressure field interpolation routine (\texttt{getprs.f90}) computes pressure fields from model variables using hydrostatic and non-hydrostatic formulations. This is critical for vertical interpolation and coordinate transformations.

\textbf{Hydrostatic Pressure Calculation:}
\begin{equation}
\frac{\partial p}{\partial z} = -\rho g = -\frac{p}{R T_v} g
\end{equation}

Integrated from surface upward:
\begin{equation}
p(z) = p_{surface} \exp\left(-\int_0^z \frac{g}{R T_v(z')} dz'\right)
\end{equation}

\textbf{Non-Hydrostatic Pressure:}
For high-resolution applications:
\begin{equation}
p = p_{hydrostatic} + p_{non-hydrostatic}
\end{equation}

where the non-hydrostatic component is computed from the vertical momentum equation:
\begin{equation}
\frac{\partial w}{\partial t} + \vec{u} \cdot \nabla w = -\frac{1}{\rho}\frac{\partial p_{nh}}{\partial z} + g + F_z
\end{equation}

\textbf{Terrain-Following Coordinates:}
In sigma coordinates ($\sigma = (p - p_{top})/(p_{surface} - p_{top})$):
\begin{equation}
p(\sigma) = p_{top} + \sigma(p_{surface} - p_{top})
\end{equation}

\subsection{Sigma Level Interpolation}

The sigma level interpolation routine (\texttt{getsiga.f90}) performs vertical interpolation between different sigma coordinate systems used by various atmospheric models.

\textbf{Sigma Coordinate Definition:}
\begin{equation}
\sigma = \frac{p - p_{top}}{p_{surface} - p_{top}}
\end{equation}

\textbf{Interpolation Formula:}
\begin{equation}
\phi(\sigma) = \phi(\sigma_k) + \frac{\sigma - \sigma_k}{\sigma_{k+1} - \sigma_k}[\phi(\sigma_{k+1}) - \phi(\sigma_k)]
\end{equation}

\textbf{Logarithmic Interpolation:}
For variables with exponential vertical profiles:
\begin{equation}
\ln(\phi(\sigma)) = \ln(\phi(\sigma_k)) + \frac{\sigma - \sigma_k}{\sigma_{k+1} - \sigma_k}[\ln(\phi(\sigma_{k+1})) - \ln(\phi(\sigma_k))]
\end{equation}

\section{Implementation Considerations}

\subsection{Numerical Stability}

All interpolation and smoothing algorithms must maintain numerical stability under operational conditions:

\textbf{CFL Condition for Smoothing:}
\begin{equation}
\alpha \leq \frac{\min(\Delta x^2, \Delta y^2, \Delta z^2)}{6D}
\end{equation}

\textbf{Monotonicity Constraints:}
\begin{equation}
\phi_{min} \leq \phi_{interpolated} \leq \phi_{max}
\end{equation}

\textbf{Mass Conservation:}
For conservative variables:
\begin{equation}
\sum_{all\_cells} \phi_{new} \cdot V_{cell} = \sum_{all\_cells} \phi_{old} \cdot V_{cell}
\end{equation}

\subsection{Computational Efficiency}

The algorithms are optimized for computational efficiency through:

\begin{itemize}
    \item Vectorization of inner loops
    \item Cache-friendly memory access patterns
    \item Precomputation of interpolation weights
    \item Parallel processing with domain decomposition
    \item Look-up tables for transcendental functions
\end{itemize}

\textbf{Parallel Implementation:}
\begin{verbatim}
!$OMP PARALLEL DO PRIVATE(i,j,weights)
do j = jstart, jend
   do i = istart, iend
      call compute_weights(i, j, weights)
      phi_new(i,j) = sum(weights * phi_old(neighbors))
   end do
end do
!$OMP END PARALLEL DO
\end{verbatim}

\subsection{Quality Control}

Quality control measures include:

\begin{itemize}
    \item Range checking for interpolated values
    \item Gradient limiters to prevent spurious oscillations
    \item Conservation diagnostics
    \item Convergence monitoring for iterative algorithms
\end{itemize}

\section{Error Analysis and Validation}

\subsection{Interpolation Error Estimates}

\textbf{Truncation Error Analysis:}
For bilinear interpolation:
\begin{equation}
E = \frac{\Delta x^2}{8}\frac{\partial^2 \phi}{\partial x^2} + \frac{\Delta y^2}{8}\frac{\partial^2 \phi}{\partial y^2} + O(\Delta x^3, \Delta y^3)
\end{equation}

\textbf{Round-off Error Propagation:}
\begin{equation}
E_{round} = \epsilon_{machine} \cdot \frac{\|\phi\|}{condition\_number}
\end{equation}

\subsection{Smoothing Filter Characteristics}

\textbf{Transfer Function Analysis:}
The transfer function of a smoothing filter:
\begin{equation}
T(k) = \frac{\phi_{output}(k)}{\phi_{input}(k)}
\end{equation}

where $k$ is the wavenumber.

For a simple diffusion filter:
\begin{equation}
T(k) = e^{-D k^2 \Delta t}
\end{equation}

\textbf{Effective Resolution:}
The effective resolution after smoothing is determined by:
\begin{equation}
k_{eff} = \sqrt{\frac{-\ln(0.5)}{D \Delta t}}
\end{equation}

\section{Advanced Topics}

\subsection{Adaptive Interpolation}

Advanced implementations employ adaptive techniques:

\begin{itemize}
    \item Error-based grid refinement
    \item Variable-dependent interpolation order
    \item Machine learning-enhanced interpolation
    \item Uncertainty quantification
\end{itemize}

\subsection{Multi-Scale Interpolation}

For multi-scale atmospheric phenomena:

\begin{itemize}
    \item Wavelet-based interpolation
    \item Scale-aware smoothing
    \item Nested grid handling
    \item Scale separation techniques
\end{itemize}

\section{Summary}

The interpolation and smoothing algorithms in GSI represent a comprehensive suite of numerical methods designed specifically for atmospheric data assimilation applications. These methods balance the competing requirements of accuracy, efficiency, and physical consistency while handling the diverse challenges presented by atmospheric flow fields, observational data, and model grids.

The implementation of these algorithms reflects decades of development in numerical weather prediction and data assimilation, incorporating both fundamental mathematical principles and practical experience from operational systems. The continued evolution of these methods, driven by advances in computational capabilities and scientific understanding, ensures that GSI maintains its effectiveness for current and future atmospheric modeling applications.

Key contributions of this algorithmic suite include robust handling of complex geometries, preservation of physical constraints, computational efficiency for high-resolution applications, and flexibility to accommodate diverse atmospheric modeling systems. These capabilities position GSI as a leading platform for atmospheric data assimilation research and operations.